%--------------------------------------------------------------------------
% function displays the dispersion relation for a given value of M and the
% neutral stability curves. then plots the path on which the system evolves
% and computes the eigenvalues on this path. with this, the amplitudes are
% computed. everything is plotted.
%--------------------------------------------------------------------------

function lo_make_path(k0, N_eigs)

dispersion = 0;

close all

if nargin == 1
N_eigs = 0;
end

addpath ../
p = params;

k = linspace(0, 3, 30);

% compute the dispersion relation
if dispersion
    subplot(1,2,1);
    lo_dispersion_relation(k, p.Ma_prime, 1);
    
    subplot(1,2,2);
end

% plot the neutral stability curves
lo_neutral_stab_curves(k, 1);

if dispersion
    hold on;
    plot([k(1) k(end)], [p.Ma_prime, p.Ma_prime],'k-.');
end
   

% compute h
t = [logspace(-6, -2, 20) linspace(1.1e-2, p.Tmax, 100)];
h = -(1 + lambertw(-p.beta / (p.beta - 1) * exp(((-p.beta + p.delta * t) / (p.beta - 1))))) * (p.beta - 1);

% compute time-dependent wavenumber and Marangoni number

M = (h - 1 + p.beta) / p.beta * p.Ma_prime;
k = k0 * h;

hold on;
plot(k, M, 'k--','linewidth',1);
plot(k(1), M(1), 'k*');

% compute the crossing of the path and the neutral stab curve
[t_c, k_c, M_c] = lo_comp_time(k0, 0);
plot(k_c, M_c, 'k.','markersize',12);


% compute the eigenvalues on this path
if N_eigs > 0
    ev = zeros(N_eigs, length(t));
    
    for i = 1:length(t)
        
        [tmp1, tmp2, L] = lo_comp_eigs(k(i), M(i));
        if (i == 1)
            [ef_tmp, tmp1] = eigs(L, 1, 'SM');
        end
           
        ev(:,i) = eigs(L, N_eigs, 'SM');
    end
    
    % plot the eigenvalues
    figure;
    plot(t, ev);
    xlabel('$t$','interpreter','latex','fontsize',12);
    ylabel('$\check{\lambda}$','interpreter','latex','fontsize',12);
    
    
    % compute the amplificiation factors
    lambda = ev(1,:) ./ h.^2;
    A_asy = exp(cumtrapz(t, lambda));
    
    % compute I and its maximum
    I = p.delta * log(A_asy);
    f = @(x) interp1(t, I, x);
    t_c = fminsearch(@(x) -f(x), 0);
    fprintf('max{I} = %.4e\n', f(t_c));
    
    ef = -interp1(linspace(0, 1, length(ef_tmp(:,1))), ef_tmp(:,1), linspace(0, 1, p.N));
    c0 = ef / norm(ef, 'inf');
    
%     z = linspace(0, 1, p.N);
%     c0 = 16 * z.^2 .* (1 - z).^2;
    



    lin = stab(k0, p); %, c0);

    A = max(abs(lin.y), [], 1);
    
    sig = diff(log(A)) ./ diff(lin.x);
    
    
    % plot the top eigenvalue
    ti = linspace(0, p.Tmax, 30);
    til = logspace(log10(t(1)), log10(p.Tmax), 30);
    figure
    subplot(2,1,1);
    semilogx(lin.x(1:end-1), sig, 'k'); hold on;
    semilogx(til, interp1(t, lambda, til), 'k*')
    xlabel('$t$');
    ylabel('Top eigenvalue');
    subplot(2,1,2);
    plot(lin.x(1:end-1), sig, 'k'); hold on;
    plot(ti, interp1(t, lambda, ti), 'k*')
    xlabel('$t$');
    ylabel('Top eigenvalue');
    
    
    % plot the amplification
    figure;
    plot(lin.x, A, 'k','linewidth',1);
    hold on;
    plot(ti, interp1(t, A_asy, ti), 'k*');
    
    xlabel('$t$','interpreter','latex','fontsize',12);
    ylabel('Amplification','interpreter','latex','fontsize',12);
    l = legend('$A_\mathrm{num}$','$A_\mathrm{asy}$');
    set(l, 'interpreter','latex','location','northeast');
   
    plot([t_c, t_c], get(gca, 'ylim'),'k--');
    
    
    % plot I
    figure;
    plot(t, I, 'k');
    xlabel('t');
    ylabel('I');
    
    
    % compute the turning points
    f_asy = @(x) interp1(t, A_asy, x);
    f_num = @(x) interp1(lin.x, A, x);
    
    opts = optimset('tolfun',1e-10, 'tolX', 1e-10);
    
    t_asy_A = fminsearch(@(x) -f_asy(x), t_c, opts);
    t_num = fminsearch(@(x) -f_num(x), t_c, opts);

     
    t_asy_L = fzero(@(x) interp1(t, lambda, x), t_c, opts);
    
    fprintf('t_c\t\t t_asy (eigs)\t t_asy (A)\t t_num\n');
    fprintf('%.4e\t %.4e\t %.4e\t %.4e\n', t_c, t_asy_L, t_asy_A, t_num);
    
    
    % plot the initial condition for the concentration perturbation
    figure
    plot(linspace(0,1, p.N),  c0, 'k');
    xlabel('z');
    ylabel('c_0');
    
end